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Abstract. - We present a novel method for calculating the fundamental gap. To this end, 
reduced-density-matrix-functional theory is generalized to fractional particle number. For each 
fixed particle number, M, the total energy is minimized with respect to the natural orbitals 
and their occupation numbers. This leads to a function, E^t , whose derivative with respect to 
the particle number has a discontinuity identical to the gap. In contrast to density functional 
theory, the energy minimum is generally not a stationary point of the total-energy functional. 
Numerical results, presented for alkali atoms, the LiH molecule, the periodic one-dimensional 
LiH chain, and solid Ne, are in excellent agreement with CI calculations and/or experimental 
data. 



Introduction. Rcduccd-dcnsity-matrix-functional theory has recently attracted a lot 
of attention [1-8]. Functionals of the one-body reduced density matrix (1-RDM) have been 
used very successfully in the calculation of correlation energies and dissociation curves of 
small molecules [9]. Given the success of these functionals it is interesting to evaluate their 
performance in the calculation of other properties. One particularly interesting quantity in 
this context is the fundamental gap A. It is given by the difference between the ionization 
potential and the electron affinity 

A = I — A, (1) 

where 

I = E^-E? ot , (2) 
A = E&-E»+\ (3) 

Here, E^ ot is the total ground-state energy of the neutral N electron system while E^ 1 are 
the ground-state energies of a system with charge T-l. In the chemistry literature (see, e.g. 
Ref. [10]), A/2 is usually termed the absolute hardness of a chemical species. Here, we use the 
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term fundamental gap for both finite and extended systems since the physical concept defined 
by Eq. ((T|) is the same in both cases. A method to calculate ionization potentials within 
reduced-density- matrix-functional theory (RDMFT) has recently been proposed [11]. 

Within density functional theory (DFT) one can prove [12] that the fundamental gap 
is exactly given by the orbital-energy difference Ae between the lowest unoccupied and the 
highest occupied Kohn-Sham (KS) state plus a number, A xc , which amounts to the discon- 
tinuity of the exchange-correlation potential upon adding and subtracting a fractional charge 
with respect to the TV-electron system. This discontinuity is zero for LDA and GGA [12]. 
Consequently, Ae is the prediction for the gap within these approximations. This prediction, 
however, deviates strongly from the experimental values underestimating them by typically 
30-50%. Moreover, strongly correlated materials, like FeO and CoO are predicted by LSDA 
to be metals (zero gap) while, experimentally, these materials are anti-ferromagnetic insula- 
tors [13]. Within exact-exchange, one of the variants of DFT functionals, the discontinuity 
A xc differs from zero leading to an overestimation of the fundamental gap. The results are 
very close to those obtained within Hartree-Fock theory [14]. Consequently, the calculation 
of the fundamental gap within DFT is still an open problem. 

In the present article we propose a method for calculating the fundamental gap by ex- 
ploiting reduced-density-matrix-functional theory. We derive a rigorous formula for the fun- 
damental gap and give numerical results for finite and extended systems. 

The Discontinuity of (i in RDMFT. - In RDMFT the one-body reduced density matrix 

j(r,r') = N J d 3 r 2 ...d 3 r N ^*(r',r 2 ...r N )^(r,r 2 ...r N ) (4) 

is used as the fundamental variable. Here, \S r (n...T'jv) is the many-body wave function of the 
interacting A^-electron system. As was shown by Gilbert [15], one can establish a rigorous 
one-to-one correspondence between the ground-state wave function and the one-body density 
matrix. Therefore, all ground-state observables are functionals of the 1-RDM. The main 
advantage of RDMFT compared to DFT is that the kinetic energy as a functional of 7(r, r') 
is known exactly 




Hence, writing the total energy in the form 

Etot [7] =T[i\ + E cxt [7] + E H [7] + E xc [7] , (6) 

where -E GX t[7] and -Eh [7] are the usual external and Hartree energy functionals, respectively, 
leads to an exchange-correlation energy which, in contrast to DFT, does not contain any 
kinetic-energy contributions. 

A complication within RDMFT is the absence of a Kohn-Sham system: Due to the idem- 
potency of the density matrix of all non-interacting systems it is impossible to reproduce the 
density matrix of an interacting system since the latter is always non-idempotent. Never- 
theless, one can directly minimize the total energy with respect to the density matrix. In 
practice, this minimization is replaced by a minimization with respect to the natural orbitals 
ipj and occupation numbers nj which are the eigenfunctions and eigenvalues of 7 

J dV 7 (r,r')^ (r') - n m (r). (7) 

In the minimization, a set of boundary conditions for n's and ip's should be considered. For 
the natural orbitals the only condition is orthonormality. For the occupation numbers there 
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are two conditions which are necessary and sufficient for the ensemble iV-representability of 
the one-body density matrix for integer particle number [16]. They are the particle number 
conservation X)j=i n j — N, and the condition < rij < 1. The orbital orthonormality and 
the particle number conservation conditions can be implemented via the Lagrange multipliers 
[i and 6ij and the functional to be minimized reads 

F[j] = E tot [ 7 ] - H (f^ n 3 -Nj-f^ e jk (| d?r<p*{r)<p k {r) - S jL )j . (8) 

The inequality condition < rij < 1 allows for optimal sets of n's containing the border values 
zero and/or one. We refer to these states as pinned states. For these states the derivative 
SF/Srij is not equal to zero in the range [0, 1] and one of the two borders, cither zero or one, 
is accepted as the optimal value for rij. This situation is most evident for non- interacting 
particles, where E tot is linear in all n^s. The existence of an occupation number which is 
equal to one means that the particular natural orbital exists in all the determinants with 
non-zero coefficient in the full CI expansion. This situation is rather exceptional for the 
exact wavefunction of a system of interacting particles. However, it is rather the rule than the 
exception for most of the existing approximate 1-RDM functionals for systems with more than 
two electrons. These functionals have the tendency to produce occupation numbers equal to 
one for most of the core states. Consequently, the functional derivative of F with respect to 7 
does not vanish at the minimum energy. One should mention that there are functionals which 
do not produce pinned states [17]. In DFT the situation is very different because any density 
which integrates to the correct particle number is ./V-rcprcscntable. Therefore, the functional 

F[p] = E tot [p] -p(j d 3 rp(r) - N^j (9) 

always has a minimum with SF/Sp = at the solution point so that SE to t/Sp(r) = /i. 

In order to derive a formula for the fundamental gap within RDMFT we extend the 
definition of the total-energy functional .Etot [7] to systems with fractional particle number 
M. Such systems can be described as an ensemble consisting of an N and an N + 1 particle 
state for N < M < N + 1 . The resulting ensemble 1-RDM for the fractional particle number 
M = N + 7] (TV e N, < rj < 1) is given by 

7 M (r,r') = (1 - V h N (r,r') + V1 N+1 (r,r'), (10) 

and the lowest energy of the ensemble is 

E^(l-r / )E£ t +r / E£+\ (11) 

where E^ ot and E^ 1 are the ground-state energies of the N and N + 1 particle system, 
respectively. It can be shown that, in complete analogy to the case of integer N, the necessary 
and sufficient conditions for a given 7 to be decomposed as in Eq. (JTUJ) arc: J2jLi n j = M 
and < rij < 1, where rij are the eigenvalues of 7 M [18]. 

As one can see in Fig. [TJ the derivative dE^ t /dM has a discontinuity at the integer 
particle number N, which by Eqs. |T][3]) is identical to the fundamental gap 



A = dE& 



DM 



n +v 9M 



(12) 

N-ri 
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Fig. 1 - Total energy for fractional total number of particles. 

So far, the derivation followed exactly the steps known from the generalization of DFT to 
fractional particle number [12]. In the following we prove that the Lagrange multiplier, ft, in 
([5]) is identical with the chemical potential, i.e. 

KM) = (13) 

In DFT, Eq. (fl~3|) is a trivial consequence of the variational equation SF/Sp = which implies 

that SE tot /Sp = ft [19]. Consequently, Egf-Efa. = E\p M ^]-E\p M ] = J d 3 r 5E/5p(r)\ pM (p M +" (r)- 

p M (r)) = fi(M) jd 3 r (p M+r > (r) - p M (r)) = fi{M)rj. In our context of RDMFT however, this 

equivalence is not at all obvious because dF/S^ need not vanish at the minimum energy (due 

to the states pinned at the border). To prove (|13[) . we investigate the difference 



E^ - E™ = E[^] - E\t»\ = I / d\d\> ' 



Sj(r, r 1 ) 



Using the fact that 



( 7 M +"(r,r')- 7 M (r-,r-')). 

(14) 



Srij 



Sj(r, r') 



^(r)^(r'), (15) 



oj[r,r') rij-rik 



the functional derivative of ([5]) yields 

SF 6E tot 



--p5{r-r'). (17) 



Sj(r,r') Sj(r, r 
Evaluating SF/6j via the functional chain rule, (fT7|) leads to 

6Etot =nS(r-r') + T— +V / d 3 x — 6tp ^ x) + c c (18) 

6j(r,r') t-r 1 Sn 3 8^{r,r') ^ J 5ipj(x) 8j(r,r') 

At the energy minimum, SF/Sipj = and SF/5ni = for any unpinned state /. The pinned 
states, however, contribute so that 



SE tl 



Sj(r, r') 



*- r ') + Efe^(^M- (19) 

p p 
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where we have used (|15|) and the sum runs over those pinned states for which SF/5n p ^ 0. 
Equation (fT4"| therefore reduces to 

C-^t^W^ + E^x / / d 3 rrf 3 rV; U W<(r')( 7 ^(r,/)-/(r,r')). 



(20) 



If we write the occupation numbers and orbitals of the M + ?/ system as 

M+rt M , c M+n M , c /ii \ 

n-j = + on,-, = <Pj +o(fj (21) 



we can obtain the first order corrections of ^ M+r i . The changes Snj and S(fj in Eq. (f2Tj) are 
not arbitrary. They are the changes in the optimal rij and if an infinitesimal charge rj is 
added to the system. Using the orthonormality of the natural orbitals we get 



- E& = »(M) V + y: 



SF 



Sn p + nf / d 3 r «(r)^;(r) + <5^(r)^ M (r)) 



(22) 

The second term in the square brackets is zero since the norm of the natural orbitals does 
not change and the first term is zero because the sum runs only over pinned states and for 
these Slip — 0, i.e. no states get unpinned. The states that would most likely get unpinned 
are those where the "true" energy minimum (SF/5n = 0) lies at a small distance outside the 
allowed interval. However, this distance is still finite and the infinitesimal additional charge 77 
cannot move the true minimum inside the allowed interval. This completes the proof of (|13D . 
Hence, by (| 1 2[) . we can evaluate the fundamental gap from the discontinuity of the Lagrange 
multiplier n(M) 

A= lim [/i(Af + r,) - »(M - n)]. (23) 

rj— >0 

Results. - Equation (|23p is of course proven for the exact functional. It is interesting 
to test it using approximate functionals E[y]. We employed the functional of Goedecker and 
Umrigar [3] (GU) which has the same structure as the Muller functional [1,2] with the im- 
portant difference that the self-interaction terms are explicitly removed. Our implementation 
for atoms and molecules is based on a Gaussian basis set expansion of the natural orbitals 
and uses the GAMESS [20] computer code to calculate the one- and two-electron integrals. 
We minimize the total energy with respect to both the natural orbitals and the occupation 
numbers employing a conjugate gradient scheme. Since the JV-representability conditions for 
fractional particle number are identical to those for integer [18], the generalization to frac- 
tional particle numbers is straightforward. We obtained the fundamental gap for several small 
molecules using 20 to 30 natural orbitals which were expanded in a cc-PVQZ atomic Gaussian 
basis set [21]. 

Fig. [2] shows the result of the numerical calculations for the LiH molecule using the GU 
functional. There is a step near M = 4 which is sharp at the lower edge but relatively smooth 
at the upper edge. The discontinuity of n(M) is located at a value slightly higher than 
M = 4 precisely at the point where the highest fractional occupation number crosses one and 
gets pinned. These features are due to the approximate nature of the exchange-correlation 
energy. The exact functional would show the discontinuity exactly at M = 4. In order to 
extract a numerical value for the gap from the graph in Fig. [2] we used the intersection of the 
extrapolated curve /x(M) for M > 4 and a vertical line at the position of the jump. 

The results for the fundamental gap of Li, Na, and LiH calculated with the GU functional 
are given in Tableland are in very good agreement with state of the art CI calculations. They 
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Fig. 2 Fig. 3 



Fig. 2 - The chemical potential n (in Hartree) as a function of particle number M for the LiH molecule. 

Fig. 3 - The chemical potential n (in Hartree) as a function of particle number M for solid Ne. The 
value of the gap is compared with HF, various DFT calculations [34] and with experiment [28]. 



also agree very well with experimental data. Note that in the context of DFT, these values 
are exceedingly difficult to calculate because, within standard LDA/GGA-type functionals, 
the negative ions of such small systems are not even bound. 

For periodic systems the symmetry properties of the many-body wave function imply that 

j{r + R,r' + R) = 7 (r,r') (24) 

for arbitrary lattice vectors R. This property, on the other hand, implies that the cigenfunc- 
tions of 7, i.e. the natural orbitals, are Bloch functions, Lp x fe(r), where A is a band index and 
fc is a wavevector in the first Brillouin zone [29]. Hence, the spectral representation of 7 reads 

7 (r,rO = ]Tn Afc ^ fc (r> Afc (7-). (25) 
x,k 

In principle, one should now minimize the total energy with respect to the occupation 
numbers 71.il and the natural orbitals I -P x ] i {i') as described above for finite systems. However, 
with the approximate functionals currently available, we encounter a serious difficulty: For 
the Miiller functional /x(M) does not show any discontinuity for all the systems we studied 
so far. The self-interaction correction of Goedecker and Umrigar [3] seems to be essential to 
reproduce this feature. However, in terms of Bloch orbitals, the self-interaction terms go to 
zero, i.e., the GU functional reduces to the Miiller functional if Bloch orbitals are inserted. To 
properly subtract the spurious self-interaction one has to use localized orbitals [30,31], such as 
Wannicr functions. Inserting the standard transformation from Bloch to Wannier functions, 

^,fcW=E elfcH ^( r - i? )' ( 26 ) 
R 

in Eq. ([25]) . the 1-RDM can be represented as 

7 (r, r') = £ £ 9x{ R &) o*(r' - R') ^(r - R) , (27) 
x R.R 

where ui\(r) is the Wannier function referring to band A and 

g x (R-R')=Y^n xk e^ R - R 'K (28) 
k 
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Table I - The fundamental gap for some atoms and small molecules as well as the LiH chain, and 
solid Ne from RDMFT compared to other calculations and experimental values, all gaps are given 
in Hartree. " QCI from Ref. [22], b from Ref. [23], c ionization potential from [22], Electron affinity 
from [24], d CISD using the same basis set as in RDMFT, e ionization potential from Ref. [25], Electron 
affinity from [26], 1 IDA, with CRYSTAL code [27] and the same basts set, 9 GGA, with CRYSTAL 
code [27] and the same basis set, h from Ref. [28] . 





RDMFT 


Other theoretical 


Experiment 


Li 


0.18 


0.175 a 


0.175" 


Na 


0.18 


0.169 c 


0.169 6 


LiH molecule 


0.27 


0.286 d 


0.271 e 


LiH chain 


0.64 


0.500 / , 0.509 9 




Ne solid 


0.76 


0.439 / , 0.546 9 


0.797' 1 



For systems with a gap, the Wannier functions decay exponentially. Hence, we expect that the 
products u)^{r' — R)u)\{r — R) contribute very little to 7 if R ^ R' . As a first implementation 
we neglect these off-diagonal terms altogether by making the approximation 

g x (R-R')=n x 6(R-R) (29) 

which leads to 

7 (r,r') = 5>A^(r' - R) w x (r - R) . (30) 
x.R 

Restricting the search to density matrices of the form (|3"D|) , we can then go ahead and minimize 
the total energy with respect to the Wannier functions and their occupation numbers n\ using 
the GU functional. The self interaction terms, when evaluated with Wannier orbitals, do 
not vanish and we obtain reasonable results (sec below). The restricted search over density 
matrices of the form (|30|) can be viewed in yet another way: By Eq. (|28j) . the approximation 
(|29|) amounts to neglecting the k-dependence of the Bloch occupation numbers, « n\. 
This is expected to be a good approximation for insulators. For metals on the other hand, 
the approximation breaks down completely because changes, at the Fermi surface, from 
values close to one to values close to zero within the same band. 

We implemented the minimization of the 1-RDM functionals in the space of Wannier states 
using the Wannier computer code described in [32,33]. As a first test case, we considered a 
system in one dimension, namely the LiH chain. Like in the case of finite systems, fx(M) 
shows a pronounced step. The size of this step compares favorably with the LDA and GGA 
values (see Tabic [J). Clearly, there are no experimental data available for this system but, as 
always, the LDA/GGA results are expected to be smaller than the true value. 

As a first fully three-dimensional system, we performed a calculation for solid Ne. Fig. [3] 
shows the discontinuity of the chemical potential when only the occupation numbers are 
optimized. The discontinuity in Nc appears slightly above 10 which is again due to the 
approximate nature of the exchange-correlation functional. The value of the gap, extracted 
from Fig. [3] by extrapolation, is also included in Table U and compares very well with the 
experimental value. 

In conclusion, we have presented a method to calculate the fundamental gap of finite 
systems and periodic solids within rcduccd-dcnsity-matrix-functional theory First numerical 
results were obtained using a recently proposed 1-RDM functional. For all systems studied, the 
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chemical potential shows a clear discontinuity as a function of the total number of electrons 
if all self-interaction terms are removed. The extracted gap values agree better with CI 
calculations and/or the experiment than any standard DFT or Hartree-Fock calculation. 

* * * 
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